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Abstract 

We propose the new family of the exactly solvable discrete state BCS - type Hamilto¬ 
nians based on its relationship to the six-vertex model in the quasiclassical limit both in 
the rational and the trigonometric cases. We establish the relation of the BCS Hamilto¬ 
nian and its eigenfunctions to the form of the monodromy matrix in the F-basis. Using 
the Algebraic Bethe Ansatz method for the standard BCS model with equal coupling 
the expression for the general scalar product and the determinant expressions for the 
physically interesting correlation functions for the finite number of sites which can be 
used in the numerical and analytical computations are obtained. We also compare the 
correlators with the results obtained by means of the variational method. 


1. Introduction. 

At present time the exactly solvable discrete-state Bardeen, Cooper and Schrieffer (BCS) 
model for the superconductivity [1] attracts much attention in connection with the problems 
in different areas of physics such as superconductivity, nuclear physics, physics of ultrasmall 
metallic grains and color superconductivity in QCD. The exact solution of the discrete- state 
BCS model is especially important for the study of the superconducting correlations in the 
atomic nuclei and the ultrasmall mettalic grains since due to the finite number of particles the 
description in terms of the grand canonical ensemble [1] (in contrast to the microcanonical 
ensemble) is obviously not correct. For all of the above mentioned problems and also for the 
many-body problems of fermions with the long-range interaction, it is desirable to find the 
integrable BCS- type Hamiltonians with the attraction of the Cooper pairs, depending on the 
momentum (on the indices of sites in the discrete - state model) and on the occupation numbers 
on the other sites. Note that at present time the study of both the continuum and the discrete 
BCS model in the case of the equal spacing of the energy levels is not completed since the 
possibility of varying the number of pairs with the filling of the part of the energy levels by the 
single- electron states was not considered. Another important problem which motivates the 
study of different modifications of the BCS Hamiltonian is to find the realistic integrable BCS 
- type models which apart from the interaction of pairs, include the interactions describing 
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the disintegration of pairs i.e. the hopping of two single electrons (fermions) to another energy 
levels independently of whether they form a Cooper pair or not. In this case the determinant 
expressions for the correlators obtained in the present paper could be useful as well. The 
studying of the excitations of different kind in the case of the microcanonical ensemble for the 
BCS model is also an interesting problem. For the case of the continuum limit of the BCS 
model as well as for the case of the microcanonical ensemble the analytical and the numerical 
calculations for various correlation functions are important. 

The eigenvectors and the eigenvalues of the discrete BCS Hamiltonian where first con¬ 
structed by Richardson [2],[3] in the context of nuclear physics. Later the model was applied 
also to the case of the Bose gas [4]. The norm and the simplest correlation functions have 
been studied in ref. [3], [4], Recently the integrability of the BCS model was proved by 
Cambiaggio, Rivas and Saraceno [5]. The set of the commuting operators which have much 
in common with the Gaudin magnets [6] was constructed. The connection of the Gaudin 
magnets with the six- vertex model was pointed out by Sklyanin [7]. The continuum limit of 
the BCS model for the equally spaced distribution of the energy levels was first considered 
by Gaudin [8] and later by Richardson [9], who also performed the numerical calculations 
[9], [10], and was found in agreement with the original (variational) BCS treatment [1]. Re¬ 
cently the solution of the model based on the off-shell Algebraic Bethe Ansatz construction 
[11] was given by Amico, Falci and Fazio [12], The numerical calculation of the correlation 
functions for the finite systems was presented by Amico and Osterloh [13] using the method 
of calculation of the scalar products, based on the generating function, proposed by Sklyanin 
[14], Nevertheless the calculations are quite involved and restricted to the systems with the 
small number of sites N. Therefore the explicit determinant expressions for the correlation 
functions of the model are highly desirable. Recently the connection of the BCS model with 
the twisted inhomogeneous six-vertex model was elaborated in ref’s [15] and [16], where the 
possibility of computation of the correlators with the help of the Algebraic Bethe Ansatz 
method was pointed out. However, the suitable determinant formulas for the correlators have 
not been obtained. The new multiparameter families of the BCS -type models connected 
with the trigonometric six-vertex model where proposed in ref’s [17], [18] and studied in more 
detail in ref. [19]. The connection of the BCS model with the conformal held theory and the 
modified Knizhnik-Zamolodchikov equations [20] have been studied by Sierra (for example, 
see [21]). The discrete BCS Hamiltonian was generalized to the case of arbitrary degeneracy 
of different energy levels of the Cooper pairs [6] and to the case of the Dicke model [22] (for 
example, see [23]). 

The main goal of the present paper is the calculation of the correlators in the BCS - type 
models with the help of the methods developed in the context of the Algebraic Bethe Ansatz 
approach to the six-vertex model. We present the new determinant expressions for the basic 
correlation functions of the BCS model. We also review and propose the new approaches to 
the construction of various integrablc BCS - type Hamiltonians. 

In the first part of the paper we briefy reveiw some of the recent results on the BCS -type 
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models, in particular, the derivation of the generalized BCS models from the six-vertex model 
both in the rational and the trigonometric cases in a way which is similar to that of ref’s [15] 
- [18]. We present the new way of construction of the BCS - type models with the interaction 
depending on the lattice sites by means of considering the limit of the transfer matrix spectral 
parameter t —> oo. The several examples of the Hamiltonians with the position - dependent 
interaction are presented explicitly. We establish the relation of the BCS Hamiltonian and its 
eigenfunctions to the form of the monodromy matrix in the F-basis which also leads to the 
generalization of the BCS model to the case of the interaction between pairs depending on 
the occupation numbers on the other sites. Let us stress, that the formalism of the F- basis is 
used not to obtain the known formulas for the scalar products for the six-vertex model, but 
to develop the new formalism for the construction of various BCS- like integrable models. 

In the second part of the paper we present the results on the calculation of various correla¬ 
tion functions for the BCS model using the Algebraic Bethe Ansatz method for the six-vertex 
model. First, we obtain the simple expressions for the scalar products and the formfactors of 
the model taking the quasiclassical limit in the corresponding expressions for the six-vertex 
model. Then, to obtain the expressions for the correlators, we use the commutational relations 
for the operators directly in the quasiclassical limit in order to reduce the problem to the cal¬ 
culation of the scalar products. We obtain the new determinant expressions for the different 
correlation functions which are usuful for the numerical and the analytical calculation of the 
correlators. Let us stress, that our results for the correlators are different from the results ob¬ 
tained previously in ref’s [13], [16], and allow one the much more simple numerical evaluation 
of the correlators, since in all cases the correlators are represented as the determinants or a 
finite sum of the determinants. 

The content of the paper is as follows. In Section 2 we propose the new family of the 
exactly solvable discrete BCS - type Hamiltonians based on its relationship to the six-vertex 
model in the quasiclassical limit both in the rational and the trigonometric cases. We present 
the examples of the Hamiltonians with the interaction between the pairs depending on the 
energy levels. We also review the results [17] on the BCS - type integrable models with the 
double set of parameters. In Section 3 we establish the relation of the BCS Hamiltonian and 
its eigenfunctions to the form of the monodromy matrix in the F-basis which also leads to 
the generalization of the BCS model to the case of the interaction between pairs depending 
on the occupation numbers on the other sites. The expression for the general scalar product 
and the determinant expressions for the physically interesting correlation functions for the 
finite number of sites which can be used in numerical computations are obtained in Section 
4. The comparison with the expressions for the correlation functions obtained with the help 
of variation on the parameters is presented in the Appendix B. Thus we show that no special 
technique of the type proposed in ref. [14] for the calculation of correlators for the BCS model 
and for the Gaudin magnets is required. The results for the correlation functions for the 
Gaudin magnets are also interesting from the theoretical point of view. The results obtained 
can be useful for studying the correlation functions for the XXZ quantum spin chain. For 
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completeness some of the results on the diagonalization of the BCS Hamiltonian are also 
presented in the Appendix A. As an additional application of the solution of the BCS model we 
present the solution of the modified Knizhnik-Zamolodchikov equations [20] for the correlators 
of the conformal held theory (SU (2) WZW- model) in the Appendix C. We present in the 
Appendix D the brief review of the Gaudin’s solution of the BCS model in the thermodynamic 
limit. 


2. BCS Hamiltonian and the Six-Vertex model. 

The BCS Hamiltonian with an arbitrary parameters e t has the form: 

N 

H = J2^ni-gS + S~, (1) 

i= 1 

where n t = bf b, is the number of the hard-core bosons (electron pairs), S + = Xa bf, S~ = J2i h 
and the coupling constant g is positive, which corresponds to the attraction. 

To solve the Schrodinger equation for the Hamiltonian (1) and its generalizations (see 
below), instead of the usual transfer matrix of the six-vertex model with twisted boundary 
conditions, we consider the following transfer matrix for the twisted inhomogeneous six-vertex 
model 

Z(t, (0) = Tr 0 (S 10 (6, t)S 2 0 (6, t)... S N 0 (6v, tj) , (2) 

which can be equivalently represented as the trace of the monodromy matrix in the auxiliary 
space 0, 

Z(t, {£}) = Tr 0 (T„(i)), T„(i) = ( ® ( ( [j 

where & are the inhomogeneity parameters. The matrices Sio(£i,t) are equal to 

Siofat) = K 0 S i0 (Zi,t) 

where A'o is the twist matrix and 5)o is the usual S-matrix of the six-vertex model obeying 
the Yang-Baxter equation S'i2S'i3S23 = <S , 23 <S'i 3 <S'i 2 : 

e v/2gN 0 \ q 

Q e -n/2.gN J r t 2 ) = h ~ t 2 + - (^ 1 ^ 2 ) 

(we denote (cricr 2 ) = ’52 a °'i 0 '2i a = X :V^ Z ) f° r the rational case. Due to the well known 
property of the matrix R 00 / = Poo'^oo' (Pi 2 = |(1 + (cr^)) - is the permutation operator) 
[A'oPA; Poo'] — 0 (it does not matter if the standard twist angle is imaginary or real) the ma¬ 
trices Sio(£i,t) obey the usual Yang-Baxter equation R 00 iS i0 S i0 ' = S i0 iS i0 R 00 > and the transfer 
matrices (2) commute at different values of the spectral parameters [Z{t)\Z(t')} = 0. The 
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Algebraic Bethe Ansatz method (for example, see [24]) can be readily applied to the mon- 
oclromy matrix defined in (2) to obtain the spectrum of the BCS model and evidently leads 
to the same results as for the transfer matrix with the usual twisted boundary conditions 
Z(t) = Tr 0 (A' 0 T 0 (f)), where T 0 (t) is the usual monodromy matrix constructed from the ma¬ 
trices Sio(£i — t) and K 0 is the diagonal matrix corresponding to the total twist angle g/2g. 
However the construction presented above allows for another generalization which will be 
considered later. Considering the limit rj —► 0 we get up to a factor (A — t): 

Sio = 1 + {r)/2gN)<j* + ( 77 / 2 ^ - t))(a 0 <Ji) + 0{g 2 ), (3) 


and retaining the terms of order 0(g 2 ), we obtain the following Hamiltonian depending on 
the parameter t: 


H{t) 


1 ^ 1 * 1 v 1 
2 9 hi (t- &) ai + 2 - o) 


(Wj) > 


(4) 


where af can be substituted by the number of the hard-core bosons af = 2n t — 1. Since the 
operators H(t) commute at different values of the parameter t, one can obtain the set of the 
commuting operators, which generalize the Gaudin magnets [6], taking the limit t —> 


Hi = --m 
9 


1 ( a t (J 1 ) 

2 i^i (6 — 6 ) 


[H t ; Hj] = 0 [H; H,] = 0. 


(5) 


Note that since the commuting operators can be defined up to an additive constant, we choose 
the operators (5), omitting the constant term in the operator af = 2n x — 1. Note also that in 
these equations the operator (cycr.,) can be represented through the hard-core boson operators 
bf, bi as 


(crjCTj) = 2 (bfbj + bjbi) + (2 m - 1)(2 rij - 1). 

These operators also commute with the Hamiltonian (1) which was first found in ref. [5]. In 
fact, considering the limit t —> 00 in eq.(4) and retaining the terms of order 1/t 2 , we 
obtain exactly the Hamiltonian (1) with e t = ^ (in what follows we omit the normalization 
factors and the additive constants depending on the total number of the hard - core bosons). 
Alternatively one can consider the following linear combination of the operators Hi to obtain 
the same Hamiltonian: 

N N 


H = -gJ2 ZiHi = E - 9S + S~, (6) 

i =1 i =1 

which coincides with the expression (1) up to an additive constant, depending on the number of 
bosons M. The following correspondence between the spin and the hard-core boson operators 


is used: S ± = JA E, Sf = 


-rr a k = 
5 ^i 


± 


= bf~(bi) = 4(erf ± ierf). One can also consider an 


arbitrary linear combinations of the operators Hi with the coefficients e* A A- Then the 
following generalization of the BCS Hamiltonian appears: 


N 

H = E e i‘ n i 

i= 1 


9 


tT. 

Kj 


(& - 0) 


( (7idi 


(7) 
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Let us note also that using the different notations the Hamiltonian (4) can be represented in 
the form: 

1 at _ 

H = - E tiUi + E e * e i ( a i a j) ■ ( 8 ) 

9 i =1 i<cj 

Thus, in general, one can construct the Hamiltonians with the interaction between the Cooper 
pairs depending on the momentum of the pairs. Another Hamiltonian of this type can be 
obtained by considering the derivative of the transfer matrix over the spectral parameter 
t. In this way the Hamiltonian takes the form similar to (8) with an extra factor e~ 1//2 + 

_i /o 

€j 7 contained in the sum of the second term. In general, since H(t ) eq.(4) depends on 
the additional parameter t and commute at different t, [H(t);H(t')\ = 0, one can take the 
derivatives over t and consider the linear combination of the operators H^it) at arbitrary t 
(H = Z n C n HW(t)). Thus, we obtain the new Hamiltonians of the BCS type (8) with the 
infinite number of additional parameters C n . The other way to obtain new Hamiltonians is 
to consider the decomposition of H(t) in the powers of l/t n at large t. In this way the new 
models with the dependence on the additional set of parameters can be obtained. Clearly, 
the similar procedures can be applied for the case of the trigonometric (hyperbolic) six-vertex 
model. 

To obtain the eigenvalues of the BCS Hamiltonian one can start with the well known 
procedure of diagonalization of the transfer matrix (2)- the Algebraic Bethe Ansatz method 
(for example, see [24], [25]). The eigenstates are represented as 

\</>(t i )) = B(t 1 )B{t 2 )...B(t M )\0) i (9) 


where the parameters t\.. Am obey the system of Bethe Ansatz equations which do not depend 
on the spectral parameter t: 

e -r)/g TT ( U-t,*- h/ 2 \ = tt f tj-tg- rj \ ' 

0=1 \ti - Ca + h/2 / 0=1 \ti - tg + T) ) ' 



Decomposing eq.(10) to the first order in the small parameter r/, we obtain the Richardson’s 
equations 


N -| Mi i 

— 2 = --• 

0=1 U So g^i U tg g 


( 11 ) 


The corresponding eigenvalue for the Hamiltonian (1) is obtained as the limit at t —> oo of 
the eigenvalue of the transfer matrix Z(t): 


M 


E = T.t i- 

2—1 


The eigenvalues of the conserved operators Hi (5) are easily evaluated from the eigenvalue of 
the transfer matrix (for the rational case) 


A (t) = e~ v/29 [] 

a 


£ a -t- r)/2 
Ca-t 


n 


tj-t + rj 
U-t 


+ e^n 

a 


jg-t + g/2 

ia~t 


n 


t-tj + rj 
t-U 
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which in the limit 77 —^ 0 up to the terms of order r/ 2 and t —> £, produces the eigenvalues of 
the operators Hi (5): 


Et 2 5 & - u) + ^ & - 6) 

(here we take into account the definition of the constant term corresponding to the operators 
Hi (5)). In the quasiclassical limit the operators B(t), C(t ) reduce to the operators E^f) 
obeying the commutational relations which can be easily obtained from the basic commuta- 
tional relations for the elements of the monodromy matrix (Yang-Baxter equation): 



N CT ±,Z 

B(t)-* £ + (t), C(t)-> £-((), £*■’(«) = E ' 


£-(«);£+(«') 


= 2 


E 2 (f) - E 2 (t') 

/ - r : 


£*(*); £±(0 


YA(t) - S ± (f') 


(13) 


= =F 


t - t' 


Note that one could built up the eigenstates of the Hamiltonian (1) directly in terms of the 
operators S ± (t) [12] (see also the Appendix A). Thus we see that instead of the off-shell Bethe 
Ansatz approach [11] used in ref. [12] one can use the usual on-shell Bethe Ansatz equations 
for the six-vertex model. 

Let us consider the trigonometric S- matrix. Repeating the steps leading to the Hamilto¬ 
nian (4) in the rational case (i.e. first taking the limit 77 —> 0) we obtain in the trigonometric 
case the following Hamiltonian: 


1 N 1 

H (t) = -7T- 53 ct s( t - &K + 5] —77 — TV~r +—FT 
2sm F - fi)sm (t - Q 


Yij + 53 ctg (t - &)ctg(t - £j)Aij, (14) 

i<j 


where the following notations are used: 

= bfbj + b~jbi, = mrij + (1 - n*)(l - rij). 


That is the one - parameter generalization of the Hamiltonian (4). First, one can proceed in 
the following way. Taking the limit t £ t , one readily obtains the trigonometric analogs of 
the commuting operators Hi (5) 


Hi = - — Tii 

2 9 


S U n(&- 6 ) 


Y u + 


-Ai. 


tg(& - 6 ) 

and considering the linear combination e t H, , we get the Hamiltonian 

€i — e 


N 


H 2 9 g>‘ + g4&-V e+ § 


tg(?i— j4,j ’ 


(15) 


(16) 


The second possibility is to consider the limit t —> 00 in the hyperbolic version of eq.(14). 
Then we get the following integrable Hamiltonian ( 7 * = H - 1 ): 


1 N 

H = ~Y a 53 l 2 i n i + 53 

y i= 1 i^j 
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The omitted terms vanish in the limit considered. 

Using the different notations the last Hamiltonian can be represented in the form: 


H = 



b?bj +bj bi 


( 17 ) 


Note that the last expression is different from the Hamiltonian (8) obtained in the case of 
the rational six-vertex model. As in the rational case the equation (14) is more general than 
eq.(16) and taking the derivatives over t one also can obtain the various generalizations of the 
Hamiltonian (17). For the hyperbolic case which is obtained by considering as an imaginary 
parameters the Hamiltonian (16) in the limit g —» oo can be considered as the Hamiltonian 
of the open spin chain with the long-range interaction. Let us mention also that the other 
possibility to construct the exactly - solvable BCS- like Hamiltonians is to consider the various 
bilinear combinations of the form H = J2i,j fijHiHj with different coefficients f t j. 

Let us comment on the form of the eigenstates and the equations for the eigenvalues for 
the trigonometric case. The analogs of the operators S ± (t) in this case are 

s± (f) = E ■ JT a t, 

i sm u - 6 ) 

and the analogs of Richardson’s equations in the trigonometric case, which can be obtained 
from the trigonometric version of the equations (10), are 


N 

E 


=! tg (U ~ £ a ) 


M 

2 V—- 


t a ) 


1 

5 
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which represent the conditions for the common eigenvectors for the operators (15) and the 
Hamiltonian depending on the double set of parameters (17). Considering the quasiclassical 
limit of the eigenvalue A(£), we obtain, omitting the factor n«sin(£ Q — t), the expression 


A(() = n Sill(t, .r * U 1 + e 


f=i sin (U - t) 


rj/2g TT sin(f } ; - t - Tj) JX- sill(^ - t + //) 

id sin (ti - t ) ^ sin(^ Q - t) 


where the definition of the parameter t differs from the definition in the rational case by the 
shift t — > t + 77/2 and considering the limit t — > we get 

Ei = l/2g — ^2 ctg(b) — Ci)- 


Here the difference with the expressions (12) is due to the different definitions of the operators 
Hi (note that = Ej + E = (1/2)(1 + (cqcrj))). 

Let us mention also the generalization of the construction proposed above. Consider the 
monodromy matrix constructed from the S- matrices S i0 (£i,t) = S i0 (£i,t) where the twist 

matrix now depends on the site i: 

e v c i/2gN q \ 

0 e ~VCi/2gN I ’ 


rv(b _ 

l\ 0 ~ 



where c t are an arbitrary parameters. Clearly in this case the Algebraic Bethe Ansatz method 
can be used in the usual way. Let us denote Co = YS <k- Then, hrst taking the limit rj —»• 0, 
we will obtain the same Hamiltonian with the coupling constant depending on the constant 
Co only. However, considering this limit in the case of Co = 1 and the parameters q are of 
order q ~ N/rj, then, to obtain the Hamiltonian of the type (1), one should commute all 
the matrices A’//' 1 to the left, which effectively leads to the gauge-like transformation for the 
operators bf, bi, which enter the operators 1 S' ± in the Hamiltonian (1). 

The discrete BCS model (1) can be generalized to the case of arbitrary degeneracy of 
different energy levels corresponding to the energies of the Cooper pairs which is equivalent to 
the case of an arbitrary spin Sj assigned to the site i with the energy q. Previously, in Section 
2 , the case s* = 1/2 was considered. One can show that the limit s* —» oo for the number of 
sites €i corresponds to the special generalization of the Dicke model. Thus, using the general 
six-vertex model, the eigenfunctions and the eigenvalues of the Hamiltonian can be obtained. 
Clearly, instead of the elementary S'-matrix one can take the following Lax operator in the 
equation (2): 

Lio(£i -t) = £i-t + r) (a Si ), 

where Sf, a = x,y,z are the spin operators with the value of the spin Sj, (Si) 2 = Sj(s* + 1). 
Considering the quasiclassical limit of the transfer matrix (2) we obtain instead of (1) the 
Hamiltonian 

N 

H = £>S?-sS + S-, ( 18 ) 

i -1 

where S ± = J2iL\Sf. Clearly, in terms of the initial BCS model the integer parameters 2sj 
correspond to the degeneracies 2s* of i-th level with the energy e t . The construction of the 
eigenstates is the same as in Section 2 and the equations (11) become 


N 


E 

a=l 


2 Si 


M i 

2 E^— 



(19) 


while the eigenvalues of the Hamiltonian have the same form E = U- The eigenstates 

can be constructed also with the help of the Gaudin operators E^t) = J2 i= 1 with ^ = €*. 

For the BCS model with the equal degeneracies 2s at each site the equations (19) take the 
form 


N 

( 2 »)£ 


M 


CK=1 


Co 


2 E 


a^i 


' tj — t r 


1 

9 


( 20 ) 


Considering the limit s —> oo of the spin at the single site, and using the expressions of the 
spin operators at this site through the Holstein- Primakoff bosons, [</>; <p + \ = 1, S z = (j) + 4> — s, 
S + = 4> + (2s - 0 + 0) 1/2 , S~ = (2s - (p + (j)) 1/2 (j) , we get after rescaling the parameters the Dicke 
Hamiltonian [22]: 

N 

H = uj(j) + (j) + 6^ - g ( S + 4> + S~4> + ) . (21) 

i =1 
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Thus the spectrum and the eigenstates of this Hamiltonian can also be obtained using the 
appropriate limiting procedure of that of the transfer matrix of the six-vertex model. Re¬ 
markably, this model can be generalized to the case of several species of the oscillators pf, pi 
with different frequencies ay. 


3. Connection with the monodromy matrix in the F- basis. 


Let us show how the above results can be obtained using the operator expression of the 
monodromy matrix in the F -basis, the basis obtained with the help of the factorizing operator 
F introduced in ref. [26]. One can construct the operator F = T\ jv which diagonalizes the 
operator Apt) [26],[27],[28] ( A F {t ) = F~ 1 Apt)F). Using the notations 


c{t) 


pit + rj) ’ 


b{t) 


(p(v) 

pit + 7]) ’ 


where <p{t) = t for the rational case and <p{t) = sin(t) for the trigonometric case, the diagonal 
operator A F (t) has the following form: 


N 

A F (t ) = ]Q (c(& - t)(l - Hi) + Hi) . 

i= 1 


( 22 ) 


Let us briefly mention some of the properties of the operator F. The explicit form of the 
operator F is 


F 12 ...N = FiF 2 ■ ■ ■ F n , Fi = (1 - hi) + Tihi, 


(23) 


where hi is the operator of the number of particles (spin up) at the site % and the operator T n 
is given by the equation 


T n ^n+l,n^n+2,n 


.s 


Nn- 


One can obtain the following formulas for the matrix elements of the operator F [28] in the 
following form: 

F{m}{n} = ({ra}|£(fm)£(fn 2 ) ■ ■ • 5 (£n M )l°)> 

where the sets of coordinates { m} and {n} label the positions of the occupied sites. The similar 
expression can be obtained for the inverse operator F -1 . Apart from diagonalizing the operator 
Apt), the operator F is the factorizing operator [26] in the following sense. For any permutation 
of indices cr G Sn ( Sn - is the group of permutations) we have the equation F = F a R a , where 
F\ 2 ..n = F a i a 2 ..aN (including the permutation of the inhomogeneity parameters £,;) and R° N is 
the operator constructed from the S- matrices defined in such a way that for the permutation 
of the monodromy matrix = To, 0 - 10 - 2 ..o-tv we have = (i? CT ) _ 1 T 0 i? fT . For the particular 
permutation cr({n}) such that ol — n\,... aM = um (tu < n 2 < ■ ■ ■ < tim) the factorization 
condition is represented as F(T ' 7 ^ n }))- 1 — T ni ..T nM . To prove the factorizing property of the 
operator (23) it is sufficient to consider only one particular permutation, for example, the 
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permutation (i,i + 1 ), since all the other can be obtained as a superposition of these ones for 
different i. One can show, that F = which evidently proves the factorization 

property. 

The matrix elements of the operators B(t) and C(t ) in the F - basis: B F (t) = F~ 1 B(t)F 
(and the same for C(t )) have the following form 

B F (t) = ^ cr\ b(£i - t ) Yl (c(£ k - t)(c(£ k - &)) _1 (1 - n k ) + n k ) . (24) 

i k^i 

C F (t) = Y, &(& - t) n (c(& -<)(!- n k) + (c(& - . (25) 

i k^i 

The operators (24) and (25) are quasilocal i.e. they describe flipping of the spin on a single 
site with the amplitude depending on the positions of the up-spins on the other sites of the 
chain. It is easily seen that in the limit 77 —^ 0 these operators reduce to the operators E ± (f) 
(13). The operator D r (t) can be found, for example, using the quantum determinant relation 
and has a (quasi)bilocal form. In fact, the following operator identity for the elements of the 
monodromy matrix (2) can be derived: 

D(t)A(t - 77) - B(t)C(t - 77) = n ((t - + 77 / 2 )(* - 6* - 77/2)), 

a 

which is readily transformed to the F- basis. From this relation the explicit form of the 
operator D F (t) can be obtained. However perhaps the simplest way to obtain D F (t) is to use 
the well known basic commutational relation following from the Yang-Baxter equation: 

|B((); C(q)] = lUi~T) - D(q)A(t)). 

Considering this equation in the F- basis in the limit q —» 00 , we obtain the following expression 
for the operator D F {t ): 

D F (t) =A F (t)+ [B F (t)]C~] , 

where the operators A F , B 1 * are given in (22), (24) and the operator C~ = linig_ + 00 (g/ 77 )C' F (g) 
equals: 

c~=~Yi b > n (( x - n k) +(c(& - ^)) _i n fc ). 

i k^i 

To hnd the physically interesting BCS-type models it is not necessary to consider the quasi- 
classical limit. For example, considering the limit t —» 00 for Z b (t) we obtain (omitting an 
overall factor rj/t and an additive constant depending on M) the Hamiltonian of the form 

H = [B + -C~] , 

where the operator B + equals 

B + = - £ bt n (wfc + {c{ik - 6)) _1 ( 1 - n k j) , 

i k^i 
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and the operator C~ was defined above. This Hamiltonian has a (quasi)bilocal form and, 
apart from the set of the parameters {£}, depends on the additional parameter rj. 

So the transfer matrix in the F-basis Z F (t) = A F (t ) + D F (t) represents the Hamiltonian 
of the BCS -type model with the varying interaction of Cooper pairs (depending on the 
occupation numbers on the other sites) even without taking the quasiclassical limit. It is 
easily seen that in the limit rj —>• 0 the results of the previous section are reproduced in such 
a way that the second term in eq.(4) corresponds to the term D h (t). Note also that since 
the pseudovacuum state is invariant with respect to the action of the factorizing operator 
F|0) = |0), the eigenfunctions of Z F {t ) have the form 11; B F (t;)|0). As it was mentioned 
above, in the quasiclassical limit the operators B F (t) and C F (t ) eqs.(24), (25) reduce to the 
operators X + (t) and X _ (t) defined in eq.(13). Concluding this section, let us stress once more 
that the Hamiltonian Z F (t), generalized by including the twist angle, contains the terms of 
the form Xu QTi; and the terms of the bilocal form ~ bf bj with the amplitude depending on 
the occupation numbers on the other sites, and reduces to the BCS Hamiltonian (1) in the 
quasiclassical limit with the corresponding twist angle. 


4. Correlation functions. 


Here we derive the analytical expressions for the simplest physically interesting correlation 
functions: (0|rq|0), (0|6d”|0), (0|S' + S'~|0) and (0|(<7j<x,)|0). For simplicity we restrict ourselves 
to the rational case although the similar formulas can be easily obtained for the general 
trigonometric (hyperbolic) models. First, we consider the following scalar product: 


F m ({A}, {*}) = (0\C(X 1 )C(X 2 )...C(X M )B(t 1 )B(t 2 )...B{t M )\0), (26) 


where {A} and {t} are the two sets of parameters, the set {£} satisfies the Bethe Ansatz 
equations and {A} is an arbitrary set of parameters. According to the connection with the 
six-vertex model revealed in Section 2, to obtain the expression for the scalar product in the 
case when the set of the parameters {t} satisfies the equations (11) one can use as a first step 
the known formula for the six - vertex model in the case when the parameters {t} satisfy the 
usual Bethe Ansatz equations (10) [28], [29], [30], [31] (see also the direct proof in ref. [27] 
and in a different way in ref. [32]), and then decompose it in powers of i] (extract the leading 
power rj 2M ). The formulas for the XXX- spin chain are 


Sm ({A}, {t}) 


1 

hi i<jv-i ~~ tj) rij<i(A; 


— det(Mij(t, A)), 
A j) *3 


where the matrix My equals: 


A) 


V f a (Aj)/ + (Aj) _ / (Aj) \ 

0 U ~ A j) V U ~ A j + 77 U- A j - 77 ) 


(27) 


(28) 
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where the following notations are used: 

M 

/ ± ( a ) = n (*<* - x± v) * 

a=l 


and a( A) = n« c(£ a ~ A). Note that in eq.(26) the operators B(t), C(t ) are normalized in 
such a way that their quasiclassical limit is given by the operators (13). From the equation 
(28) taking the limit A* —» t t one can easily obtain the formula for the norm of the Bethe 
eigenvector: 

N M (t) = t 3 + n) det [tv ( t )] ; 

lh ftiU-tj) v 

where the matrix N tJ can be represented in the form: 




2 rj 


(tij + v){Uj ~ V) ’ 


^ 3), 


N u = 


d_ 

dt, 


In (a(ti)))-^2 — 

a^i WO 


2 rj 


+ V) (tai ~ V) ’ 


where t l3 — ti — tj and a(t) = e 2ri ^ J fla c(£a — t), a(t ) —> 1 + 2 rj/g + rj J2 a (t ~ £<*) 1 at rj —» 0. 
Extracting from the last expression the term of order r ) 2M , we get the norm of the eigenstate: 


(0(f) 10(f)) = det [-Wij(f)], 

*7 


(here (0(f)10(t)) = SM(t,t)) where the matrix iV^ is given by 


*« = E 

a 


l 

(ti - ■Ca ) 2 


E 



% 



* 7 ^ 3- 


(29) 


The formula for the norm was first derived by Richardson [3]. Note that the norm does not 
depend on e* explicitly. Similarly the general scalar product (26) is given by the formula (27) 
with the matrix of the following form: 


Mij(\,t) = 


(U - A j) 2 


n(f<* 


a Saj k^i 


1 1 

( A j ffc) 9 


(30) 


Thus the scalar product (26) with one Bethe eigenstate is obtained in the form given by the 
equations (27), (30). Equivalently, since the product in the matrix (30) depends only on the 
number of the column j this scalar product can be also represented in the form 


Sm({ A}, {f}) 


FI i,j(ti X j) 

n i<j(ti — tj) 


YT d.et (M tj (A, t)), 

x j) 13 


where the matrix M l3 {\, t) is given by 


Mij(\t) 


1 

(f* - x i) 2 


E 


i 

(A, - £a) 


2E 


1 

(A j — tk ) 
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However for the calculation of the correlators in order to consider the limit A * —> fj for some 
of A i, it is more convenient to use the formulas (27), (30). The formula for the norm can 
be obtained directly by taking the limit A, —> U in this expression for Sm(A, f) for the BCS 
model. Since the set {t} satisfies the Richardson Bethe Ansatz equations (11) we obtain from 
(30) in this limit the matrix elements 


m u -»• n 

OL^l 



u p) ’ Mij —► n (t 

a^i j aft 




*7 ^ 3- 


Again, since the product which depend only on the number of the column j can be written 
in front of determinant, taking into account the products in front of the determinant (27), we 
obtain exactly the norm (29). For completeness, let us present another equivalent expression 
for the scalar product with one Bethe eigenstate which can also be used for the derivation of 
the norm and can be useful for taking different limits in the process of the calculation: 


Sm( A, t) — 


rii^; (ti A j ) 


El i<j(ti tj) rij<i(Aj A j) ij 


———d.et , 


where the new matrix M i3 is equal to 


Mij(\,t) 


0 tj Aj) ( ST' 

(U ~ Aj) 2 Aj - 6 


2E 




5 


where the hrst term in the numerator also depends only on the index j. 

Now one can calculate the physically interesting correlation functions using the Algebraic 
Bethe Ansatz method. Although the expressions for some of the correlators can be obtained 
by the method of variations over the parameters [3], [6], which we present in the Appendix 
B, the determinant expressions obtained with the help of the Algebraic Bethe Ansatz method 
are different in the form (although equivalent to the variational ones). At the same time the 
direct Bethe Ansatz method allows for the computations of some of the correlation functions 
(such as ( bfbj ), for example) that are not accessible by the variational method and can be 
useful for the computation of the correlators in the different BCS- like models. 

Let us proceed with the calculation of the simplest correlation function (rij) along the 
lines of ref. [27] and using the technique developed above. First, we consider the action of the 
operator A(£j) at the state [0(A)) (9) where A* are not necessarily satisfy the Bethe Ansatz 
equations. Due to the symmetry of the problem, one can consider the site £i. One can use 
the formulas for the scalar product with Bethe eigenstate (27), (28) taking subsequently the 
quasiclassical limit. Note that due to the symmetry of the Hamiltonian (1) here it is not 
necessary to use the general formulae of the Quantum Inverse Scattering method for the six- 
vertex model [33]. However, the most simple way is to represent the eigenstate directly in 
terms of the operators (13), 


|0(Q) = E + (f 1 )S + (f 2 )...E + (f M )|0> 


(31) 


14 



and use the formula for the scalar product (30). To calculate the expectation value (0(t)|ni |0(t)) 
we act by the operator n\ on the state (31) using the following relation: 


niS + (t) = S + (t)ni + 


(t ~ 6) 


bj, m|o) = o. 


The result can be represented in the form which allows one to apply the expression (30) for 
the scalar product: 

ni\<f>(t)) = Y, [f .\ i) ^(C-6) (s + (ti)...S + (C) ...S + (t M )) | 0 ), 

where the operator S + (C) replace the operator £ + (C) in eq.(31). The factor (£ — £i) —> 0 in 
this formula is cancelled by the corresponding term in the denominator in the first sum in 
the expression (30) for the matrix M l} at A* = (. Next, one can use the following theorem. 
Consider the determinant of the sum of the two matrices det jj (+ a^-), where the second 
matrix a>ij = Cj0j - is the matrix of rank 1. Then we have: 


M 


det (Mij + dij) = det(My) + cy. det(M^), 
ij U /, | '■> 


(32) 


(k) 

where the matrices M ,h differ from the initial matrix only by the substitution of its k-t\i 
column by <f>f 

' — (1 — 5jk)Mij + Sjkfpi- 

Applying this theorem to the case of the average (0(A)|ni|0(t)), the determinant in eq.(30) 
transforms into the sum of the determinants in (32) with 


,. = m (fc) = 1 

* (**-6) 2 


Thus, we get the expression 
(0(A) |ni | (/>(£)) = 


rn^-6), c fc = n 

a a^k 


Arv A/ 


A Q — £i / (A k - 6) 


rw*i tj) n_/<i (Aj Aj) 


(det (.!/,, + Hij) - det (Mi , 


where the matrix H tJ = Cj4>i in the limit A* —> t, t equals 


H ij {A tj 






and the matrix M.\j is given by the equation (30). Taking the limit A, 
expression, we finally obtain the formula 


ti in the whole 


(0(t)|ni|0(f)) = det (N i:j + H^) - det(iV ii ), 


(33) 
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where iVy(t) is the matrix of the norm (29) and the rank-one matrix Hij equals 


ij (ti-6) 2 ' 

To find the expectation value one should divide this expression by the norm of the eigenstate 
1 4>(t)) which is given by det(iVy). Thus the expectation value (n*) is represented as a ratio 
of the determinants. This formula can be used in the numerical evaluation of the occupation 
number in the case of finite system (for finite N). 

The same method can be used to obtain the determinant representation for the two- 
point correlation function ( bfbj). Note that for this correlator the final expression was not 
obtained in ref. [3]. Due to the symmetry of the problem it is sufficient to calculate the average 
(0(f)|0(f)). First, we act by the operator b\ on the state \<f>(t)) using the formula 

6 1 S + (f) = z + (t)h + ——(1 - 2 m). 

[t - SlJ 

Then the action of the operator b\ produce the state (see also the Appendix A): 

bi\m = E ( 2 + (*i) •••(*)••• s + (*»)) |o> 

where we denote by (i), (j) the absence of the operators E + (C), E + (C) in the last product. 
Clearly, the action of the operator b^bi produces the state with two terms (see the above 
formulas and the Appendix A for more details) where the first term is equal to 

bp>A<Kt)) m = E y (C - W (2 + («i)... E + «)... E+(«„)) |o>, 

where the operator S + (^) is substituted instead of i-th operator S + (C). Proceeding in the 
same way as for the average (rq) we get the similar expression for the hrst term 

(0(f)|fcjAi|0(t)) (1) = det(Nij + H\f) - det(%), (34) 

where the rank-one matrix is equal to 

M i) = 1 (tj - 6) 

y (u - h ) 2 & - ZiY 

The second term for & 2 ~&i|</>(£)) contains the double limit of the form 

—2 E 77 TTT, -FT Jim j™ (C. - 6KC2 - 6) (s+(i 1 )...S+(Ci)...S + «2)...) |0>, 
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which means that the two columns should be replaced in the resulting determinant. Namely, 
repeating the above procedure and taking into account the form of the matrix My (30) we 
come to the expression 





1 

fi) 


/ (4-6)(4-fr) 



5 


where the expression in the parenthesis comes from the factor which can be written in front 
of the determinant in eq.(30) and the matrix N^' l} equals N l3 with the exception of the two 
columns k,l, which are equal: 


N (k,l) = (4 ~ £l) pr(k,l) = (tl ~ 6) 
(U - 6) 2 ’ (4 ~ 6) 2 ’ 


In the equivalent form the last expression for the second term is 

1 

k<i 


(0(f)|^&i|0(f)) (2) = -2]T 


V(4-t z )(6-4)y v 

where the columns /e, / of the new matrix are 


det (JV^>) , 


!Cr(fc.O _ 4*0 _ (4 ~ £ 2 ) 

^ “(4^p 


N^ l) = 


,(o _ (4 - 6) 


(4-6) 2 


Note that in the particular case of the average {bfbi) due to the first term we immidiately 
get the expression for (ni), since the determinant in the second term is equal to zero clue to 
presence of two identical columns in the matrix Thus, the remarkably simple expression 

as a sum of determinants is obtained. 

Let us perform the similar calculations for the average ( riiTij ). Considering the action of 
the operator n^rti on the state j0(f)), we get the expression similar to the second term for 

Hbi\(/>(t)): 


n 2 ni| 0 (f)) = ^ - — lim lirn (Ci - 6 XC 2 - 6 ) (s + (fi)...S + (Ci)...S + (C 2 )...) |0). 

i^j (4 - 4 l )(tj ~ 42) Cl-’■Si C2-& v 2 

Repeating the calculations performed for the average ( bfbj ), we obtain the following sum of 
the determinants 


(0(f) |n 2 ni |0(f)) = X) 


'(4 -6XR-6)' 


v(4 -4)(6 -6) 

where the columns k, l does not depend on the indices k and l: 


det . 


Nit" = 


dfc) 


(4-£ i) 2 ’ 


Nt l) = 


di) 


(4 -6) 2 ’ 
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The last expression can also be represented in the equivalent form since, taking into account 
the factors in the parenthesis, one can modify the the columns of the matrix AT ’ . Namely, 
we have 

(0(f) |n 2 n! \<P(t)) = E 77-7777-FT det (E?’°) » 

(4-f/)(6-4i) v v ' 


r( fc .O 


where the two columns of the new matrix Ay. '' are 


jCr(k,l) _ 

I '^ ik 


dk) 


(4-6) 

(t.-6) 2 ’ 


EE = 


00 


(fr - 6) 

(4 -6) 2 


Note that the above expression for ( riiUj) is obtained for i ^ j. Clearly, at i — j the 
average coincides with (n0. Since both averages (bfbj) and (nirij) are found, the correlator 
((oyoj)) can also be calculated. It is straightforward but lengthy calculation to show that 
the above expressions lead to the expression for (( a^Gj )) found by the variational method and 
presented in the Appendix B. Thus, the remarkably simple determinant expressions for the 
pair correlators of the model suitable for the numerical calculations even at the sufficiently 
large N are obtained. 

As for another physically interesting correlator - the expectation value (S + S~), it could 
be obtained either using the correlator (bfbj) obtained above or directly with the help of the 
representation of the average as a certain limit of the scalar product. However in fact it is 
sufficient to calculate the average (Xy 6 n i) using the formula (33) which leads to the final 
result 

(0(t)| E& n *l0(f)) = det (Nij + Hij) - det(Ajj), (35) 

i 

where Hij is again the rank-one matrix of the form 


A = £ 

a 


6 

(0 - 6) 2 


Clearly, due to the form of the Hamiltonian (1), from the equation (35) one can obtain the 
expression for the average (S + S~). We will show below that, in fact, the sum of eq.(35) and 
the expression for ( S + S~) obtained by the different method gives the total energy E. 

Let us note that since the average should be divided by the norm of the state in all of the 
above mentioned cases the expectation values of the operators can be represented in the form 

det (Stj + (AT 1 #)..) , 

where is the Kronecker symbol and the matrix H tJ stands for one of the three matrices 
introduced above (33), (34), (35). Since in the case of the averages (n0 and ( S + S~) (or, 
equivalently, (Xa6ffi)) 4 ie rank-one matrices H y, Hij depend only on the first index i, the 
last expression can be represented in a different form with the help of the identity det(<5jj+yj) = 
1 + Xi- Namely for the first average we get 

<t> = E ( Ar_ 1 )p- h' h = (ti-^y (36) 
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From the equation (36) one readily obtains for the average ( S + S ) = (Y: X' n i) / 9 ~ E/g 


E(a} — y (tv _i 




y 


‘'31 


= E 




- ZaY 


( 37 ) 


in agreement with the expression (35) obtained with the help of the different method. From 
the equation (36) one can check that the total number of particles is YY) = M which can 
be easily seen from the following property of the matrix N t f. 

E N * = fib), f(u ) = E 77 \ \ 2 ’ 

a a S a) 


see eq.(29), or equivalently 1,; = Y YEq/(E- 

Now let us turn to the calculation of the average ( S + S~) in a direct way. With the help 
of the formula S + = lim^^ C£ + (C) we get using the equations (11) that the operator S + S~ 
acts on the state | (f>(t)) as 

S+S-\m) = -Y. lim C (E + (ii)...E + (C)(i)-S + (iM)) |0), 
g Y C^oo v J 


where S + ((7) stands on the i-th place (see equation (42) in the Appendix A). Here the factor 
1 /g appears due to the equations (11). Then taking the limit ( —» 00 we finally obtain 


(cj)(t)\S + S | (f>(t)) = det (N i:j + H i:j ) - det(Nij), 


with 

lJ 9 2 ’ 

where the second factor 1/g is due to the last term in the matrix (30). The matrix Hij does 
not depend on the indices at all. Dividing this expression by the norm and repeating the steps 
described above, we obtain the following formula for the average ( S + S ~): 

< sh - s -> = 4 e (38) 

9 Xj V 13 

in agreement with the result based on the variational method. This form is analogous to the 
expressions (36), (37), for the other correlators. With these three expressions one can check 
the consistency of the above formulas. Namely the sum of (37) and (38) should be equal to 
the total energy of the system E. To prove it we use the following equation for the matrix 
Nif 

Y N- t = V_ — _ - 

A/, 2 ’ 

a a \ L i S a) y 

which can be easily obtained using the equations (11). Multiplying both sides of this equation 
by the matrix N~ 3 , we immediately get the energy E as a sum of (37) and (38). 
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In the cases, when the expressions for the correlation functions can be obtained with the 
help of the variation over the parameters of the total energy and the integrals of motions 
[3], [6], the results coincide with that obtained above as the ratio of the determinants. That 
can be easily proved for the averages (rq), ((cqcr.,)), ( S + S ~) using the formula for the matrix 
inverse to the norm matrix N t] (29). We discuss the correspondence of this two approaches 
in more details in the Appendix B. 


Conclusion. 

In conclusion, for the applications to the realistic fermionic systems, it would be interesting 
to find the realistic discrete - state integrable Hamiltonians with the interaction of fermions 
containing the terms describing the breaking of the Cooper pairs. It is possible that the models 
proposed in the present paper can be generalized to this case. Let us mention that even in 
the framework of the BCS model (1) one can include the terms describing the interaction 
of the single- electron states in the following way. Namely, consider the lattice Hamiltonian 
constructed from the fermionic operators c+ (c iCT ) with a =j, j= 1,2: 

H = Y e i( n l* + n 2i) + Y v.ij (. SjSj ) - ()Y ( K b j + b t b i) f 

i i<j i<j 

where ni a = c^c^, bf = and Sf = ^(cf afci) (the sum over the spin indices 1,2 is 

implied). Since for this Hamiltonian the number of double- occupied sites is conserved the 
eigenstates are given by the superposition of the eigenstates of the BCS Hamiltonian and the 
eigenstates corresponding to the second term Vij(SiSj). Here the coefficients V l} are not 
necessarily the constant but for an integrable model can correspond to any of the integrable 
quantum spin chains (for example, the XXX- spin chain or the Haldane - Shastry spin chain in 
its trigonometric [34] or hyperbolic [35] versions). Note also that using the method presented 
the Hamiltonians of the discrete- state BCS- like models related to the exactly solvable t — J- 
models of different symmetry [36] and the models based on the different Lie algebras both in 
the rational and the trigonometric cases can be constructed. The main shortcoming of these 
models from the point of view of the applications to the description of electrons close to the 
Fermi- surface is the presence of the single-electron hopping terms. In this context the study 
of Gaudin magnets for the different Lie algebras can be useful. 

We have shown that for the calculation of the correlators for the Gaudin magnets and the 
BCS model no special technique for the calculation of the scalar products [14] is required. 
Using the determinant expressions for the correlation functions obtained in the present paper, 
one could hope to ford the analytical results for the correlators in the thermodynamic limit 
for the systems with different density of energy levels. The determinant expressions for the 
correlation functions can also be useful for the numerical evaluation of correlators both in 
the case of the BCS model with the fixed number of pairs and the general BCS, which takes 
into account the existence of the single - occupied energy levels for the excited states. In 
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conclusion, let us note, that the similar determinant expressions can be obtained both for 
the BCS and the Gaudin magnets models based on the different Lie algebras, which is an 
interesting problem from the theoretical point of view. 
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Appendix A. 

Here we present the most simple and beautiful procedure of diagonalization of the Hamil¬ 
tonian (see for example [6]) based on the construction directly in terms of the operators (13). 
We look for the eigenstate in the form: 

|0(f)) = H+(f 1 )H+(f 2 )... E+(t M )|0> (39) 

We use the following commutational relations for the operator E + (t) which can be easily 
proved directly or derived from the general formulas (13): 

niE + (f) = E + (t)m + 

&iE + (f) = E + (f)6 1 + —(1 - 2m). (40) 

First, using the first relation, we obtain the formula for the action of the operator n a to the 
state (39) |0(t)) = 11): 

n a \t) = £ 7t 4-T E+ ( ( ‘) •••(*)••• 2 + («*,)y |0), 

i S a) 

where the sign (i) means the absence of the operator E + (tj) in the product. Then considering 
the sum ]T a £ a n a we obtain the following expression for the first part of the Hamiltonian: 

J2^n a \t) = -S + ]T|t(i)) + 1 1), (41) 

OL i Vi/ 

where we denote by the sign \t(z)) the state (39) without the single operator E + (fj). Next, 
consider the action of the term — gS + S~ on the state (39) using the second commutational 
relation (40). First, considering the action of the operator b \, we obtain the formulas: 

i'll t) = E (E + («,)... (1 - 2m )... £+(«„)) |0) = 

? (i^6) 1 *»> - 2 g fe-Jfe-f.) ( E+(tl) ■'' (i) " w) ''' E+( ‘" ) ) ^ |0> - 
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After the simple algebraic transformations, using the definition of the operator £ + (f), we 
obtain the following expression: 


s-\t) = E (E TT^TT) l*W) + 2 E 7—^-7 ( 1 * 0 )) - l*(0». 

i \ a fb Za) / i<j Z j) 


which finally leads to the result: 


s + s-\t) = E | E rzr - 2 E rwr 1 s+ li(i)). 

i \ a SO b r a 


(42) 


Thus combining the equations (41) and (42) we get 


H\m = mt ))-sE (E T~rr~ 2 E—X + I S+ I<W), 

i \ a Sa Q ^j b fa y 


and obtain the eigenvalue E = Xa U and the condition of the cancelation of the “unwanted” 
terms S + \t(i)) - the equations (11): 


E 



2E 

a^i 




l 

9 


It is straightforward to find also the eigenvalues of the conserved operators (5) using this 
method. In fact, performing the similar calculations, we finally obtain the formula for the 
action of the operators H, L (5) to the state | (f>(t)) (39). For the operator 


H i = 


1 

9 


n i 


;E 

H i 


(°i°j) 

(6-60 


we get the expression 


Hi\t) 


ly. 1 I y- 1 

2 a (6 - fa) + i (ti ~ 6) 


10- 


E 



2E 


i 


ti tfc 



1 

E - 6) 




The condition of cancelation of the “unwanted” terms bi\t(i)) is equivalent to the equations 
(11) and the eigenvalues of the operators Hi are given exactly by the equation (12), however, 
the last expression for H\ t) is valid for the state of the form (39) for an arbitrary sets of the 
parameters {t} and {£}. 

Note that the similar expression for H{\t) can be easily obtained also for the trigonometric 
case even without the detailed calculations. In fact, it is clear, that the analog of the last 
formula has the form 


H 1 \t) = E 1 {t,£)\t) + 'Ef( t i) 

i 


l 

sin(f* - fi) 


b i\t (*)>, 
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where Ei(t,£) is the eigenvalue of the operator Hi and the function /(U) equals 

f(ti) = Y Ct §(^ - fa) - 2 Y Ct §(^ - ta) + l/g, 

a a 

so that the trigonometric Richardson’s equations have the form f(ti) = 0. The coefficient in 
front of the second term in the formula for H\\t) can be fixed from the known expression for 
the action of the operator ri\ contained in H\. 


Appendix B. 

Here we compare the results for the correlators obtained with the help of the Algebraic 
Bethe Ansatz method with the results obtained by means of the variation over the parameters 
[3], [6]. Variation over the parameters £* and g gives the equation: 

A SU 5g 

a (U - &*) 2 (*i - ^) 2 0 (ti - ^) 2 9 2 ' 

For the case of the averages (rq), ( S + S ~) i.e. varying over the parameters g this equation 
leads to 

N iJ St i = Y faptitp + Sg/g 2 , 

3 g 

where N tJ is the norm matrix (29) and the matrix (f ) t 0 = l/(t t — ^) 2 . The solution of this 
system of equations gives the variation of the energy 5E = Yi &U which determines the above 
averages. One can easily see that for the correlators (rq), ( S + S~) and also (Xafdh) the 
expressions obtained in Section 4 starting from the determinant expressions are reproduced. 

In order to calculate the average ((cqcrj)), we consider the variations of the eigenvalues 
(12) of the commuting operators (5) over the parameters 

m 1 i 

SZj 2 (£; - Q* {atCrj) ' 

We get for the eigenvalues 

5Ei _ 1 1 ^ 1 5ti 

Hi ~ 2 (6 - Q 2 v 


where according to the formulas of Section 4 the matrix 5ti/5£@ equals: 

Stf i\ 1 




«(o - W 


Using these formulas we easily obtain for the average the expression 

1 n 1 


((Wj)) = 1 - 2 (& ~ £j) Y 


l,k 


(ti - w 2 


W 1 ) 


lk (tk-tj) 2 ' 


23 



which can be shown to be in agreement with the determinant expression given in ref. [ 6 ]. In 
fact it is straightforward to represent this average as a ratio of the determinants: 

((cr^)) = 1-2 (f* - Q 2 det (Rij)/ det(iVjj), 

where the matrix is (M + 1) x (M + 1)- matrix with the indices i,j = 0,1,... M and with 
the matrix elements equal to i? 0 o = 0 and 

Rkl = Nkl, RkO — 77 -yrx, Roi — 77-777 

(4 - Cj) (*i ~ Si) 

for k, l — 1... M. Thus the determinant expression is obtained. 


Appendix C 

Using the formulae obtained above, we present here the solution of the modified Knizhnik- 
Zamolodchikov equations for the vector- valued function [$(£)) of N variables £i.. .£jv corre¬ 
sponding to the SL( 2 ) algebra: 


d 1 , 1 v- ( cr i cr i) 

7 36 g n ' 2^J(& -6) 


!*(«)> = o. 


(43) 


The modification corresponds to the term ~ 77 and at g equal to infinity and the additional 
parameter 7 = 1 we obtain the usual Knizhnik-Zamolodchikov equation [20] for the correlators 
of the conformal held theory (WZW- model). We consider the rational case although the same 
procedure can be performed for the trigonometric case. Our presentation follows the solution 
[11] and make use of the eigenstates of the BCS- model (or twisted six-vertex model in the 
quasiclassical limit). However, we present here the solution which is different (in the form) 
from the approach [11], and based on the off-shell Bethe Ansatz equations for the model (1), 
or, equivalently, on the equations for the six-vertex model after the quasiclassical limit was 
already performed. For simplicity consider the case 7 = 1. 

Following ref. [11] let us seek for the solution of the equation (43) in the form: 

[$(£)) = j>dt x(t,0 10 (f)), 

where | (f>(t)) (39) is the eigenstate of the Hamiltonian (1), considered as a vector- valued 
function of the variables t t and and x(£,£) is some function to be specified below. We 
denote by the symbol j> dt the integration over the variables t\.. Am in the complex plane 
over some closed contours C\... Cm- The particular solution of the equation depends on the 
choice of the contours. Taking the derivative of this vector- valued function over the variable 
£1 gives 

= fdt 1 4>(t))+ jdt (44) 
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Differentiating the eigenstate | </>(t)) we obtain 


where the state | t(i)) is the state (39) without the single operator E + (fj). Since we assume the 
closed contours of integration over 0 and the integration by parts can be applied, the second 
term in the equation (44) can be represented as 

To rewrite this term we make use of the obvious similarity between this formula and the 
expression for given by the last formula in the Appendix A, which is valid for an 

arbitrary parameters 0 and In fact, one can represent this equation in the equivalent form 
as 

£/(«*)773tAI««} = - ffi) 10. 

i (A si) 

where 

/w=e—^- 21 :— 

a L i Sa L i L k y 

If one can choose the function y(t,£) in such a way that (d/d tj)x(t, £) = f(ti)x(t,£) and the 
last term takes the form 

dt X(t,0 (Ei(t,£) - H±) |0(t)), 

where Ei(t,£) is the eigenvalue of the operator H\ (12), and the operator Hi is obviously does 
not depends on 0, then we immediately find that the function |$(£)) is the solution of the 
equation 

(A + ^)wo> = °. 

provided the function x(t,£) satisfies the equation (d/d£i)x(t,£) = Ei(t,£)x(t,£)- Thus we 
get the following system of the equations for the function y(f, £): 


d_ 

d£i 


■xiut) 



AS 

a^i 





d_ 

dti 


x{t,£) 



2E 


1 


t'i tk 



x(t,0- 


One can see that this two sets of the equations are compatible and the function %(f, £) is given 
by 

x (t, o= e^/s n& - tj )- 2 n (e« - ^) i/2 mu - u 

i<j a<f3 i,a 
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Remarkably, the variation of this function over the variables U gives the equations (11) 
if one assumes that U are at the stationary points corresponding to the function x(f,£) 
{{d/dti)x(t,£,) = 0). For the parameter 7 ^ 1 the similar expression for the function x(f,£) 
can be obtained (the powers of the factors in the pruducts should be divided by 7). In the 
similar way the solution of the equations based on the different Lie algebras can be obtained. 
In the trigonometric case the solutions of the SU (2) KZ - equations have the similar form: 

r\ r\ 

0 = Ei(t, Z)x{t, £), 0 = f{U, t, £), 

where the expressions for E^t, £) and f{U,t, £) should be substituted by their trigonometric 
analogs presented above: 


Ei = 1/2 g - X]ctg {tj -&), f(U ) = (U - £a) - 2 ^ ct s(^ - t a ) + 1/g. 

7 a a 


From these equations the explicit solution for the function xi/,0 i n the trigonometric case 
can be obtained. 

Let us mention that in the rational case, except the commuting differential operators 
d/dt/i — Hi, in the SU{ 2) case there is an extra commuting differential operator g 2 d/dg — H, 
where H is the Richardson Hamiltonian (1), (6) and g is the corresponding coupling constant: 



— H 


d_ 

d£i 


-H 


0 


(for example, see [37], [38]). In conclusion, let us mention that the integration contours 
Ci,... Cm, which define the solution of the modified KZ - equations, are not necessarily the 
closed contours, which do not intersect the branch cuts of the integrand in the complex plane. 
The only condition is that the integral over the total derivative of the form 


E 


dt 


d_ 

dt. 


x(t,0 


{U - £*) 


KW)) 


0 


for each a — 1,... N, which follows from the derivation presented above. 


Appendix D 

Here we present the analytical solution of the BCS model in the continuum limit (i.e. 
in the limit N —» 00) for the case of the equal- spacing distribution of the energy levels £, : , 
or for the density of energy levels p{£) equal to unity at some interval which can be chosen 
as £ G (—1,1). Although the solution of the equations (11) for the continuum limit has 
been considered previously [6], [9], with the help of the electrostatic analogy, and the final 
solution in agreement with the result of the variational BCS treatment (1) was obtained, some 
questions remained obscure. For instance, the choise of the ansatz for the complex electric 
held with a branch cut along the line is not well understood. Thus the additional arguments 


26 



and approachs for this problem are highly desirable. First of all, note that the equations (11) 
can be represented as the conditions of the minimum of the “energy” functional d/dti&(t) = 0, 
where the roots ti are considered as a positions of charges at the two- dimensional complex 
plane interacting through the Coulomb potential and subjected to the homogeneous electric 
field of the strength — 1/g, directed along the real axis, 

$(£) = ln lb - 6 *| “ 2 In| U -tj\ + -J2 Ret*. 

i,a i<j 9 i 

The potential between the like charges t; of the value +1 corresponds to the repulsion, while 
the their interaction with the points £; with the charges — 1/2 is repulsive. 

Let us denote by h(z) the holomorphic (which depends only on the coordinate z at the part 
of the complex plane without the charges) complex electric field defined in such a way that 
the integral of h(z) over some closed contour C in the cornlcx plane is equal $ c dzh(z) = 2-Jiiq, 
where q is the total electric charge enclosed by the contour C. Clearly, the unit charge, 
centered at the origin, gives the complex electric field equal to \jz. 

First, let us suppose that in the continuum limit the roots ti form some curve T, symmetric 
with respect to the real axis, with the continuous complex density of roots along the curve 
R(t) (we assume here that all roots t, appears in a complex pairs). Let us denote by a and 
b = a* the endpoints of the curve T, and by C- the contour enclosing the curve T. We also 
assume that the curve T does not intersect the domain of the charges £ at the real axis. The 
discontinuity of h(z) at V is equal to the density R(t), A h(z) = 2niR(z), which by the Cauchy 
theorem means that h(z) is the analytic continuation of R(z) from the curve T. In particular 
that means that the total number of particles is § c dz h(z) = 2niM and the total energy is 
§ c dz zh(z) = 2iriE 

The Gaudin’s assumption is that the field h(z) has the branch cut along the line (a, b) of 
the form 

tx(z) = \J(z — a) (z — b) J ^770(0, 

where the function </>(£) can be fixed from the condition that the residues of the field at the 
points £ should be equal to — 1 / 2 : 

^K) = (l/2)(K-a)K-6))- I/2 . 


Since at infinity 


oo the field should be equal to —1/g, we have the equation 


/ 


d£0(£) = 


1 

V 


The conservation of the total number of particles gives the equation f r dtR(t ) = M, and the 
total energy equals E = f r dttR(t). 

Alternatively, one can write down the equations (11) in the continuum limit in the following 
form: 




t-t 


— 2 dz 


m 

t — z 


1 

9 
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where the integal over the curve T can be transformed to the integral over the contour C as 



Substituting the ansatz for h(z), we obtain 


1 

9 






2 / <%</>(£) 


dz a)(z - b ) 

c 2m (t — z)(£ — z) 


and transforming the integration contour C into the two contours, the first one is around the 
domain of £, and the second one is at the circle at the infinity, we get 






This equation leads to the two equations, the first one is the equation for the function 
and the second one is for the integral / d/</(/) = 1/2 g. Applying the same transformation of 
the integration contour to the integrals M = § c (dz/2m)h(z) and E = § c (dz/2ni)zh(z), and 
substituting the value a(b) = ±iA, we obtain the gap equation, the equation for the particle 
number and the energy in agreement with the BCS theory. To invetsigate the 1/N corrections 
Richardson [9] has derived the closed Riccati type integro-differential equation for the electric 


held 


11 11 
h(z) = V -V- 

i Z-U 2 « Z~£a g 


in the continuum limit, using the same operation with the contour integration in the complex 
plane as above. However, at present time, the solution of this equation (without using the 
ansatz for the held h(z)) is absent. The form of the curve T can be obtained from the condition 
that the component of the electric held along the curve should be equal to zero for the point 
at the curve T. One can imagine the curve T as a metallic plate of the special form with the 
endpoints a, b. At each point near this plate the vector of the electric held has a direction 
perpendicular to the plate. In particular that means that the electic held at the points a, b 
is equal to zero, which is fulfilled for the ansatz for h(z). In other words, the curve T can be 
found from the condition that it should be the equipotential curve for the electric held h(z). 
One should stress that the form of the branch cut between the points a and b can be chosen in 
an arbitrary way, and for the ansatz for h(z) it is assumed that the branch cut coincides with 
the curve T. Note that one can calculate the density of charges \R(t)\ along the curve and 
obtain the form ((t ~ a)(t- b)y/ 2 f(t), where /(f) is some smooth function, characteristic 
for the matrix models. Since one can imagine the conformal mapping of T onto the line (a, b), 
which reduce the problem to the solution of the matrix model with some potential, which 
should exhibit the density of states of the same form, this can be considered as a justification 
of the ansatz. 
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